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We report calculations of the electronic structure of silicene and the stability of its weakly buckled 
honeycomb lattice in an external electric field oriented perpendicular to the monolayer of Si atoms. 
The electric field produces a tunable band gap in the Dirac-type electronic spectrum, the gap being 
suppressed by a factor of about eight by the high polarizability of the system. At low electric fields, 
the interplay between this tunable band gap, which is specific to electrons on a honeycomb lattice, 
and the Kane-Mele spin-orbit coupling induces a transition from a topological to a band insulator, 
whereas at much higher electric fields silicene becomes a semimetal. 

PACS numbers: 73.22.Pr, 63.22.Rc, 61.48.Gh 



I. INTRODUCTION 

Two-dimensional (2D) carbon crystals arc hosts for 
Dirac-type electrons, whose unusual properties have been 
studied extensively in graphene monolayers produced by 
mechanical exfoliation from graphitejii^ A close relative 
of graphene, a 2D honeycomb lattice of Si atoms called 
silicene^ does not occur in nature, but nanoribbons of sil- 
icene have been synthesized on metal surfaces.—"— Due to 
the similarity of the lattice structures, the band structure 
of silicene resembles that of graphene, featuring Dirac- 
type electron dispersion in the vicinity of the corners 
of its hexagonal Brillouin zone (BZ)ji Moreover, silicene 
has been shown theoretically to be metastable as a free- 
standing 2D crystal)^ implying that it is possible to trans- 
fer silicene onto an insulating substrate and gate it elec- 
trically. In this work wc predict the properties of this 2D 
crystal. 

The similarity between graphene and silicene arises 
from the fact that C and Si belong to the same group in 
the periodic table of elements. However, Si has a larger 
ionic radius, which promotes sp^ hybridization, whereas 
sp^ hybridization is energetically more favorable in C. 
As a result, in a 2D layer of Si atoms, the bonding is 
formed by mixed sp^ and sp^^ hybridization. Hence sil- 
icene is slightly buckled, with one of the two sublatticcs 
of the honeycomb lattice being displaced vertically with 
respect to the other, as shown in Fig. [TJ Such buckling 
creates new possibilities for manipulating the dispersion 
of electrons in silicene and opening an electrically con- 
trolled sublattice-asymmetry band gap4 In this article 
we report density functional theory (DFT) calculations 
of the band gap A for Dirac-type electrons in silicene 
opened by a perpendicular electric field using a combi- 
nation of top and bottom gates. We show that A can 
reach tens of meV before the 2D crystal transforms into a 
semimetal and then, at still higher fields, loses structural 
stability. We also determine the weak electric field at 
which electrons in silicene experience a transition from a 
topological insulator regime^iiiS caused by the Kane-Mele 
spin-orbit (SO) couplingii for electrons on a honeycomb 
lattice into a conventional band insulator regime. 

The rest of this article is arranged as follows. In Sec. 




FIG. 1: (Color online) Atomic structure of silicene, together 
with a sketch of the charge density for the highest occupied 
valence band in the vicinity of the K point. 



|TT]we report our results for the structural and electronic 
properties of free-standing silicene, and compare them 
with other theoretical and experimental results in the 
literature. In Sec. IIIII we analyze the effects of a trans- 
verse electric field on the structural and electronic prop- 
erties of silicene, and in Sec. IIVI we discuss the effects of 
SO coupling on the electronic structure, arguing that a 
crossover from topological insulating behavior to band in- 
sulating behavior must take place as the transverse field 
increases in strength. In Sec.|V]we give the technical de- 
tails of our computational methodology and demonstrate 
the convergence of our results with respect to simulation 
parameters. Finally, we draw our conclusions in Sec. IVII 



II. STRUCTURAL AND ELECTRONIC 
PARAMETERS OF FREE-STANDING SILICENE 

A. Comparison with theoretical and experimental 
results in the literature 

The lattice constant and the z (out-of-plane) co- 
ordinates of the Si atoms lying on the 2D honey- 
comb lattice were both fully relaxed using DFT (i) 
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TABLE I: Silicene structural and electronic parameters: lat- 
tice constant a, sublattice buckling Az (the diflerence between 
the z coordinates of the A and B sublattices) , cohesive energy 
-Ec, and Fermi velocity v. The calculated cohesive energy of 
silicene includes the DFT-PBE zero-point energy, which we 
found to be 0.10 eV per atom. The theoretical results are for 
free-standing silicene; the experimental results are for silicene 
nanoribbons on Ag substrates. 
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FIG. 2: (Color online) DFT-PEE phonon dispersion curves 
for silicene in zero external field and at — 0.51 VA~^. In 
both cases the calculations were performed using the method 
of finite displacements, with the atomic displacements being 
0.0423 A, in a supercell consisting of 3 x 3 primitive cells with 
a 20 X 20 k-point grid in the primitive cell. 



in the local density approximation (LDA), (ii) with 
the Perdew-Burke-Ernzerhof (PBE) exchange-correlation 
functionalfi^ and (iii) with the screened Heyd-Scuseria- 
Ernzerhof 06 (HSE06) hybrid functional^ii Our DFT 
calculations were performed using the CASTE P^^'^^ and 
VASnil plane-wave-basis codes, using ultrasoft pseu- 
dopotcntials and the projector-augmcnted-wavc (PAW) 
method, respectively. The z-coordinates of the two Si 
atoms in the unit cell (the A and B sublattices) differ 
by a finite distance Az. Our results for a free silicene 
monolayer are shown in Table HI The metastable lattice 
that we find is the same as the "low-buckled" structure 
found by Cahangirov et al^ 

The experimental results for the lattice parameter de- 
pend on the choice of substrate on which the silicene 
is growuj^i^ The extent to which theoretical results ob- 
tained for free-standing silicene are applicable to the sil- 
icene samples that have been produced to date is there- 
fore unclear. 



are summarized in Fig. [2l This convinces us that, as a 
metastable 2D crystal, silicene can be transferred onto an 
insulating substrate, where its electronic properties can 
be studied and manipulated as suggested below. 



C. Electronic band structure 

The calculated band structure of a "free" silicene layer 
is shown in Fig. [31 As expected, it resembles the band 
structure of graphene; in particular it shows the lin- 
ear Dirac-type dispersion of electrons near the K points, 
where we find the Fermi level in undoped silicene. The 
Fermi velocity v of electrons in silicene is lower than that 
in graphene (see Table H]) . Although the lattice parame- 
ters and sublattice buckling found in the different DFT 
calculations are in good agreement, our results for the 
Fermi velocity are very much smaller than the Fermi ve- 
locity reported in Ref. H. 



B. Stability of free-standing silicene 

The cohesive energy of bulk Si (including a correc- 
tion for the zero-point energy) has been calculated within 
DFT-LDA as 5.34 eVJ^ Comparing this with our DFT- 
LDA cohesive energy of silicene reported in Table U shows 
that bulk Si is substantially (0.22 eV per atom) more sta- 
ble than silicene, implying that silicene would not grow 
naturally as a layered bulk crystal like graphite. How- 
ever, by calculating the DFT phonon dispersion it has 
been verified both here and in Ref. H that the structure 
is dynamically stable: no imaginary frequencies appear 
anywhere in the BZ. The results of such an analysis 



III. APPLICATION OF A TRANSVERSE 
ELECTRIC FIELD 

A. Breaking the sublattice symmetry 

To exploit the weak buckling of silicene, we consider 
its behavior in an external electric field applied in the 
z-direction, as shown in Fig. [1] The main effect of such 
an electric field is to break the symmetry between the A 
and B sublattices of silicene's honeycomb structure and 
hence to open a gap A in the band structure at the hexag- 
onal BZ points K and K'. In the framework of a simple 
nearest-neighbor tight-binding model, this manifests it- 
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FIG. 3: (Color online) DFT-PBE band structures for silicene 
in a cell of length Lz — 26.5 A with a plane- wave cutoff energy 
of 816 eV and a 53 x 53 k-point grid; (a) in zero external 
electric field, (b) with = 0.26 VA"\ and (c) with = 
0.51 VA~^ (shown both with and without the relaxation of 
the atomic coordinates in the electric field). The zero of the 
external potential is in the center of the silicene layer. The 
dashed line shows the Fermi energy in each case and the insets 
show the spectrum near the Fermi level in the vicinity of the 
K point. 



self in the form of an energy correction to the on-site 
energies that is positive for sublattice A and negative for 
B. This difference in on-site energies A = £a — £b leads 
to a spectrum with a gap for electrons in the vicinity of 
the corners of the BZ: £± = ±^(A/2)2 + \vp\^, where 
p is the electron "valley" momentum relative to the BZ 
corner. Opening a gap in graphcnc by these means would 
be impossible because the A and B sublattices lie in the 
same plane. 



B. First-order perturbation theory 

A nai've estimate of the electric-field-induced gap in 
silicene can be made using first-order perturbation theory 
by diagonalizing a 2 x 2 Hamiltonian matrix at p — ^ 0, 



(V-kI^IV-k) (V'kI^I'/'k) 



(1) 



Here, tjj^ are the degenerate lowest unoccupied and 
highest occupied Kohn-Sham orbitals at the K point at 
Ez = 0, and z = corresponds to the mid-plane of the 
buckled lattice. This suggests a band gap which opens 
linearly with the electric field at a rate dA/dEz = 0.554 
and 0.573 eA for the wave functions found using the 
LDA and PBE functionals, respectively. 



Self-consistent DFT calculations in the presence 
of the field 



The estimate given in Sec. 1111 Bl is in fact only an up- 
per limit for the rate at which the band gap opens, since 
it neglects screening by the polarization of the A and 
B sublattices. In order to obtain an accurate value of 
the rate at which a band gap can be opened with an 
electric field, we have performed fully self-consistent cal- 
culations of the DFT band structure in the presence of 
an electric field. A typical result of such a calculation 
is shown in Fig. EJb). At small electric fields, relaxing 
the structure in the presence of the field does not have 
a significant effect on the band gap, but the screening of 
the electric potential by the sublattice polarization of the 
electron states makes a substantial difference. The DFT- 
calculated gaps are gathered in Fig. 2] The variation of 
the band gap A at K with electric field E^ is almost per- 
fectly linear for fields up to i?z « 1 VA~^. The results for 
the rate dA/dEz at which a gap is opened are shown in 
the table inset in Fig. 01 The eightfold difference between 
the self-consistent and the unscreened values of dA/dEz 
indicates that the system exhibits a strong sublattice po- 
larizability. 

Our value for the rate at which the band gap opens 
within DFT-PBE is 0.0742 eA. This is substantially lower 
than the result obtained by Ni et al.^ which is 0.157 cA. 
Part of the reason for the discrepancy is that we extrap- 
olated our results to infinite box length, whereas Ni et 
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FIG. 4: (Color online) DFT gap against applied electric field 
Ez for silicene with a plane-wave cutoff energy of 816 eV and 
a 53 X 53 k-point grid. Unless otherwise stated, the PBE 
functional was used. The box length in the z direction was 
varied from — 13.35 A to 26.46 A. The results have been 
extrapolated to the limit — > oo of infinite box length (solid 
lines) as described in Sec. IV D 41 Unscreened band gaps cal- 
culated using perturbation theory are also shown. The inset 
table shows the calculated rate at which the band gap opens. 




FIG. 5: (Color online) DFT-PBE and DFT-LDA band struc- 
tures with and without SO coupling taken into account. The 
inset shows the bands around the K-point, revealing a small 
band gap induced by SO coupling. The width of the bottom 
panel corresponds to 1/200 of the FK line. 



IV. SO COUPLING IN SILICENE 



A. SO-induced gap 



al. used a fixed amount of vacuum between the periodic 
images of the layers. Another possible reason for the dif- 
ference is that we used a plane- wave basis set, whereas Ni 
et al. used a localized basis set. An incomplete localized 
basis set would tend to undermine the extent to which 
the electrons can adjust to screen the electric field. 



We have also performed a study of the effects of SO 
coupling (which is more pronounced in Si than in C) on 
the band structure. The SO coupling term is explicitly in- 
cluded in the Hamiltonian in the DFT calculations. The 
results obtained with the LDA and PBE functionals are 
shown in Fig. [5l Both functionals predict an SO gap of 
the order of a few mcV at the K point, while the rest of 
the band structure barely differs from the nonrelativistic 
case. Our calculated LDA and PBE SO gaps are 1.4 meV 
and 1.5 meV, respectively, in agreement with the recent 
literature.— 



D. Stability of the silicene lattice in an electric field 

The narrow-gap silicene band structure shown in Fig. 
|3] persists over a broad range of electric fields Ez ■ How- 
ever, for electric fields of more than ~ 0.5 VA^^, the 
band gap starts to close due to an overlap of the conduc- 
tion band at F and the valence band at K, and silicene 
becomes a semimetal, as shown in Fig. ^c). According 
to our calculations, the buckled honeycomb crystal is still 
mctastable at this electric field, as can be seen in Fig. [21 
The main effects of the electric field on the phonon dis- 
persion curve are (i) to lift some degeneracies at K and 
M and (ii) to soften one of the acoustic branches, but 
without making the frequency imaginary. Under much 
higher electric fields, the honeycomb structure of silicene 
becomes unstable. We found that E^ > 2.6 VA~^ causes 
the lattice parameter to increase without bound when 
the structure is relaxed. 



B. Crossover from topological to band insulating 
behavior 

In the theory of Dirac electrons on the honeycomb lat- 
tice, the SO gap is accounted for by the Kane-Mele term 
describing, e.g., intrinsic SO coupling in graphenej^i The 
Kane-Mele SO coupling and the electric-field induced A- 
B sublattice asymmetry for electrons in the vicinity of 
the BZ corners K± = (±47r/(3a), 0) in silicene can be 
incorporated in the Hamiltonian 

Hk± = Vp ■ (T + AsoSzC^z + ^S.^zCrz, (2) 

where ^ = ±1 distinguishes between the two valleys, K+ 
and K_, in silicene's spectrum. Here, the Pauli matrices 
a^, (Jy, and CTz act in the space of the electrons' ampli- 
tudes on orbitals attributed to the A and B sublatticcs. 
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{ipA, tpn) for the valley at K+ and {ipB, — V'a) for the val- 
ley at K_. In Eq. ([2]), is the electron spin operator nor- 
mal to the silicene plane, and Aso and A2 are the DFT- 
calculated SO-coupling and electric-field induced gaps. 

The Hamiltonian of Eq. ([2|) generically describes the 
transition between the 2D topological and band-gap in- 
sulators. Its spectrum. 



£t± = ±^^(Aso + CA.f + «V, 

£l± = ±\l\ (Aso-^A.f + i;V, (3) 

includes two gapped branches, one with a larger gap 
I Aso + Az| and another with a smaller gap |Ago — Az|. 
At a critical external electric field E'^ w 20 mVA~^, 
Aso = A^, and the smaller gap closes, marking a tran- 
sition from a topological insulator—"— at Ago > A^ to 
a simple band insulator at Ago < A^. The difference 
between these two states of silicene is that the topolog- 
ical insulator state supports a gapless spectrum of edge 
states for the electrons, in contrast to a simple insulator, 
where the existence of gapless edge states is not protected 
by topology. However, one may expect something rem- 
iniscent of the topological properties of Dirac electrons 
to show up even in the band insulator state of silicene: 
an interface between two differently gated regions, with 
electric fields E^ and —E^ (where E^ » E^), should 
support a one-dimensional gapless band with an almost 
linear dispersion of electrons 

V. COMPUTATIONAL DETAILS 
A. Cohesive energy 

All our plane-wave DFT total energies were corrected 
for finite-basis error— and it was verified that the residual 
dependence of the total energy on the plane- wave cutoff 
energy is negligible. We used ultrasoft pseudopotentials 
throughout, except where otherwise stated. The silicene 
system was made artificially periodic in the z-direction 
(normal to the silicene layer) in our calculations. The 
atomic structure was obtained by relaxing the lattice pa- 
rameter and atom positions within DFT, subject to the 
symmetry constraints and at fixed box length in the z 
direction. The cohesive energy was then evaluated using 
this optimized structure. 

The energy of an isolated Si atom (needed when eval- 
uating the cohesive energy) was obtained in a cubic box 
of side-length L subject to periodic boundary conditions. 
We extrapolated the energy of the isolated atom to the 
limit of infinite box size by fitting 

E{L) = E{oo) + cL-^ (4) 

to the DFT energies E{L) obtained in a range of box 
sizes, where E{oo) and c were parameters determined by 
fitting. Equation ^ gave a very good fit to our data. 



We have also calculated the DFT zero-point correction 
to the energy of silicene. This is expected to be largely 
independent of the exchange-correlation functional used. 
Indeed, our calculations show that the zero-point correc- 
tion is 0.103 eV within the LDA and 0.101 cV with the 
PBE functionalii^ We used the PBE result in our final 
calculations of the cohesive energy reported in Table ID 



B. Evaluation of the Fermi velocity 

To evaluate the Fermi velocity shown in Table U we 
evaluated the DFT band structure using a 53 x 53 k-point 
grid and a plane- wave cutoff energy of 816 eV in a cell 
of length = 26.46 A. We then fitted Eq. (17) of Ref. 
[23I to the highest occupied and lowest unoccupied bands 
within a circular region around the K point; the Fermi 
velocity is one of the fitting parameters. The radius of 
the circular region was 6% of the length of the reciprocal 
lattice vectors; we verified that the Fermi velocity was 
converged with respect to this radius. 



C. Geometry optimization and phonon dispersion 
curves 

The phonon dispersion curves shown in Sec. Ill B] were 
calculated using the method of finite displacements, with 
atom displacements of 0.042 A, in a supercell consisting 
of 3 X 3 primitive cells with a 20 x 20 k-point grid in the 
primitive cell. In the results with the external electric 
field, the box length was — 19.05 A and the plane- 
wave cutoff energy was 435 eV. In the results without the 
field, the box length was ~ 13.35 A and the plane- 
wave cutoff was 816 eV. This choice was made because 
the error due to a finite box length is potentially much 
larger in the presence of a transverse electric field. 

The geometry optimization and band-structure calcu- 
lations at zero external field were performed with both 
the CASTEBi^ii^ and VASpii codes, to verify that the re- 
sults are in good agreement. This check was necessary be- 
cause it was only possible to perform the electric-field cal- 
culations with CASTEP, while for the SO calculations we 
had to use VASP. In principle the only difference between 
the calculations performed using the two codes arises 
from the Si pseudopotentials used. The PAW method^ 
was used in the VASP calculations, whereas ultrasoft pseu- 
dopotentials were used in the CASTEP calculations. As 
can be seen in Tabic [H the geometries predicted by the 
two codes agree well. We have also verified that the band 
structures are in good agreement. Finally, in Fig. [5] we 
show that the phonon dispersions obtained with the two 
codes are virtually identical when the same parameters 
are used. 

Figure [S] also demonstrates that our phonon dispersion 
curves are converged with respect to supercell size. 
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FIG. 6: (Color online) Phonon dispersion curves for silicene 
obtained with CASTEP and VASP using different exchange- 
correlation functionals and supercell sizes. The results for 
a 3 X 3 supercell were obtained with a box length Lz = 13.35 
A, a 20 X 20 k-point grid in the primitive cell, and a plane- 
wave cutoff energy of 816 eV. The matrix of force constants 
was evaluated using the method of finite displacements, with 
the displacements being 0.042 A. The results for a 7 x 7 super- 
cell were obtained with a box length L, = 15.0 A, a 12 x 12 
k-point grid in the supercell, and a plane-wave cutoff energy 
of 500 eV. The matrix of force constants was evaluated us- 
ing the method of finite displacements, with displacements of 
0.09 A. 



D. Band gap in the presence of an external electric 

field 



1. Plane-wave cutoff energy 

The convergence of the calculated band gap with re- 
spect to the plane-wave cutoff energy for a particular 
applied field is shown in Fig. [T] The gap converges ex- 
tremely rapidly. 
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FIG. 7: (Color online) DFT-PBE gap against plane-wave 
(PW) cutoff energy for silicene subject to an electric field 
of 0.257 V/A in a cell of length = 13.35 A with a 15 x 15 
k-point grid including K. 
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2. k-point sampling 



FIG. 8: (Color online) DFT-PBE gap against the reciprocal of 
the number A^^ of k-points for silicene subject to an electric 
field of 0.257 V/A in a cell of length = 13.35 A with a 
plane-wave cutoff energy of 816 eV. 



The convergence of the calculated band gap at BZ 
point K with respect to the k-point grid used in the 
self-consistent field calculations is shown in Fig. |8l The 
finite-sampling error falls off as the reciprocal of the total 
number of k points. The prefactor of the finite-sampling 
error is vastly greater when K or K' is included in the 
grid of k-points for the self-consistent field calculations. 



3. Choice of pseudopotential 

The dependence of the calculated gap on the exchange- 
correlation functional and pseudopotential is shown in 
Fig. [5] The difference between the results obtained with 



different pscudopotcntials is much smaller than the gap, 
but is not wholly negligible. The on-the-fly ultrasoft 
pseudopotential is believed to be more accurate than the 
norm-conserving pseudopotential?^ and hence we have 
used the former in our final calculations. 



4- Box length 

The dependence of the calculated gap on the length 
of the simulation box is shown in Fig. (TU] It is clear 
that this is potentially a large source of error. However, 
the L^-depcndence is reasonably well approximated by 
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FIG. 9: (Color online) DFT gap with different exchange- 
correlation functionals (LDA and PBE) and pseudopotentials 
[on-the-fly ultrasoft (USP) and norm conserving (NC)i^] for 
silicene subject to an electric field of 0.257 V/A in a cell of 
length Lz — 13.35 A with a 15 x 15 k-point grid and a plane- 
wave cutoff energy of 816 eV. 



lengths L,_ = 13.35, 18.521, and 26.459 A; the results 
were then extrapolated to infinite box length by fitting 
a quadratic in l/L^. From the magnitudes of the varia- 
tions shown in Figs. [Tl llOl we estimate the uncertainty in 
our final results for the rate dA/dE^ at which the band 
gap opens when an electric field is applied to be less than 
about 0.01 eA. 



E. Unscreened estimate of the band gap in the 
presence of an external electric field 

To evaluate the field-induced band gap using perturba- 
tion theory we used norm-conserving pseudopotentials.— 
We used a 39 x 39 k-point mesh including the K point, 
and a cell length of Lz = 26.46 A. It was verified that the 
perturbation-theory-induced rate of gap opening dA /dEz 
was converged to within 0.00002 eA with respect to k- 
point mesh and Lz- The finite-basis error in dA/dEz 
was found to fall off approximately exponentially with 
respect to the plane-wave cutoff energy, and hence we 
extrapolated our results to basis-set completeness. 
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FIG. 10: (Color online) DFT-PBE gap against the reciprocal 
of the box-length Lz for silicene subject to an electric field 
of 0.257 V/A with a 15 x 15 k-point grid and a plane-wave 
cutoff energy of 816 eV. The solid line is a fitted quadratic in 
l/Lz. 



F. Band structure with SO coupling 

The SO calculations were performed with a plane- wave 
cutoff of 500 eV and a 24 x 24 k-point grid. We checked 
that the length of the simulation box has negligible influ- 
ence on the SO gap: the gap is the same with simulation 
box lengths of 15 A and 30 A up to numerical accuracy. 



VI. CONCLUSIONS 

In summary, we have shown that a 2D layer of Si 
atoms — silicene — is a versatile material in which a band 
gap can be tuned (in a broad range of tens of meV) us- 
ing a transverse electric field Ez, while silicene remains 
metastable. At the low field Ez ~ 20 mVA~^, we ex- 
pect silicene to undergo a transition between a topologi- 
cal and a simple band insulator, whereas at much higher 
field Ez ~ 0.5 VA~^ it will undergo a transition from a 
band insulator into a semimetal. 



a quadratic in 1/Lz, allowing the DFT gaps to be ex- 
trapolated to infinite cell size if results at three or more 
different cell sizes are available. 



5. Estimates of uncertainty in our final results 

Our final results for the field-induced gap were ob- 
tained using 53 X 53 k points, not including K or K', 
a plane-wave cutoff energy of E'cut = 816 eV, and box 
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